---
title: "SPPQ-23-0075.R2"
author: "Replication Code"
format: docx
embed-resources: true
code-background: true
---

```{r preamble, echo = F, warning = F, message = F, include = F}
# Clear environment
rm(list=ls())

# Load packages
library(foreign)
library(sandwich)
library(MASS)
library(tidyverse)
library(tidymodels)
library(stargazer); library(knitr)
library(kableExtra)
library(haven)
library(xtable)
library(lme4)
library(corrplot)
library(ggmap); library(sf)
library(ggpubr); 
library(gridExtra)
library(Matching)
library(cspp)
library(urbnmapr)
library(modelsummary)
library(fixest)
options(scipen=7)
options(digits=4)
set.seed(20120630)

# Create clusterMod function for standard errors
clusterMod<-function(model, cluster)
{
  require(multiwayvcov)
  require(lmtest)
  vcovCL<-cluster.vcov(model, cluster)

  coef<-coeftest(model, vcovCL)
  #w<-waldtest(model, vcov = vcovCL, test = "F")
  get_confint<-function(model, vcovCL){
    t<-qt(.975, model$df.residual)
    ct<-coeftest(model, vcovCL)
    cse<-sqrt(diag(vcovCL))
    est<-cbind(ct[,1], cse,ct[,1]-t*ct[,2], ct[,1]+t*ct[,2])
    colnames(est)<-c("Estimate", "Clustered SE","LowerCI","UpperCI")
    return(est)
  }
  ci<-round(get_confint(model, vcovCL),4)
  return(list(coef, ci))
}

# Load in dataset for analysis
data <- read.csv("Women Friendliness Dataset.csv")
uniqcan <- read.csv("unique women candidates.csv")
```

```{r data_wrangling, echo = F, warning = F, message = F, include = F}
# Split into Dem, GOP subsamples
datdem <- data %>% dplyr::filter(CandidateDemocrat == 1) 
datgop <- data %>% dplyr::filter(CandidateDemocrat == 0)

# Split by Office Category and Party
govdem <- data %>% dplyr::filter(OfficeCat == "Governor" & CandidateDemocrat == 1)
cabdem <- data %>% dplyr::filter(OfficeCat == "Cabinet" & CandidateDemocrat == 1)
othdem <- data %>% dplyr::filter(OfficeCat == "Other" & CandidateDemocrat == 1)
govgop <- data %>% dplyr::filter(OfficeCat == "Governor" & CandidateDemocrat == 0)
cabgop <- data %>% dplyr::filter(OfficeCat == "Cabinet" & CandidateDemocrat == 0)
othgop <- data %>% dplyr::filter(OfficeCat == "Other" & CandidateDemocrat == 0)

# Split into White, Non-White
whitedat <- data %>% filter(white == 1)
nonwhdat <- data %>% filter(white == 0)

# Create race-and-party subsets
dat.demwh <- data %>% filter(white == 1 & CandidateDemocrat == 1)
dat.demnw <- data %>% filter(white == 0 & CandidateDemocrat == 1)
dat.gopwh <- data %>% filter(white == 1 & CandidateDemocrat == 0)
dat.gopnw <- data %>% filter(white == 0 & CandidateDemocrat == 0)


# To create scatterplots of votepct and wfd by party
dat.wfd <- data %>%
  dplyr::group_by(wfd_ten) %>%
  dplyr::summarise(all.mean = mean(na.omit(votepct)),
            all.sd = sd(na.omit(votepct)),
            all.obs = n())
dat.gop.wfd <- data %>% filter(CandidateDemocrat == 0) %>%
  group_by(wfd_ten) %>%
  dplyr::summarise(gop.mean = mean(na.omit(votepct)),
            gop.sd = sd(na.omit(votepct)),
            gop.obs = n())
dat.dem.wfd <- data %>% filter(CandidateDemocrat == 1) %>%
  group_by(wfd_ten) %>%
  dplyr::summarise(dem.mean = mean(na.omit(votepct)),
            dem.sd = sd(na.omit(votepct)),
            dem.obs = n())

dat.office.dem <- data %>% filter(CandidateDemocrat == 1) %>%
  group_by(wfd_ten, OfficeCat) %>%
  dplyr::summarise(dem.mean = mean(na.omit(votepct)),
                   dem.sd = sd(na.omit(votepct)),
                   dem.obs = n())
dat.office.gop <- data %>% filter(CandidateDemocrat == 0) %>%
  group_by(wfd_ten, OfficeCat) %>%
  dplyr::summarise(gop.mean = mean(na.omit(votepct)),
                   gop.sd = sd(na.omit(votepct)),
                   gop.obs = n())

dat.race.dem <- data %>% filter(CandidateDemocrat == 1) %>%
  group_by(wfd_ten, white) %>%
  dplyr::summarise(dem.mean = mean(na.omit(votepct)),
                   dem.sd = sd(na.omit(votepct)),
                   dem.obs = n()) %>%
  filter(!is.na(white)) %>%
  mutate(white = ifelse(white == 1, "White", "Nonwhite"))

dat.race.gop <- data %>% filter(CandidateDemocrat == 0) %>%
  group_by(wfd_ten, white) %>%
  dplyr::summarise(gop.mean = mean(na.omit(votepct)),
                   gop.sd = sd(na.omit(votepct)),
                   gop.obs = n()) %>%
  filter(!is.na(white)) %>%
  mutate(white = ifelse(white == 1, "White", "Nonwhite"))

# Safe, Competitive, Dangerous state comparisons
dat.dwoc.safe <- data %>%
  filter(avg_votes >= 0.55 & CandidateDemocrat == 1 & white == 0)
dat.dwoc.danger <- data %>%
  filter(avg_votes <= 0.45 & CandidateDemocrat == 1 & white == 0)
dat.dwoc.comp <- data %>%
  filter(avg_votes > 0.45 & avg_votes < 0.55 & CandidateDemocrat == 1 & white == 0)

dat.dww.safe <- data %>%
  filter(avg_votes >= 0.55 & CandidateDemocrat == 1 & white == 1)
dat.dww.danger <- data %>%
  filter(avg_votes <= 0.45 & CandidateDemocrat == 1 & white == 1)
dat.dww.comp <- data %>%
  filter(avg_votes > 0.45 & avg_votes < 0.55 & CandidateDemocrat == 1 & white == 1)

dat.rw.safe <- data %>%
  filter(avg_votes >= 0.55 & CandidateDemocrat == 0)
dat.rw.danger <- data %>%
  filter(avg_votes <= 0.45 & CandidateDemocrat == 0)
dat.rw.comp <- data %>%
  filter(avg_votes > 0.45 & avg_votes < 0.55 & CandidateDemocrat == 0)
```


```{r summary_statistics, echo = F, warning = F, message = F, include = F}
plot1.all <- ggplot(dat.wfd, aes(X = wfd_ten)) +
  geom_col(aes(x = wfd_ten, y = all.obs), size = 1, color = "chocolate4", fill = "white") +
  geom_line(aes(x = wfd_ten, y = 125*all.mean), size = 1.5, color="chocolate4", group = 1) +
  scale_y_continuous(sec.axis = sec_axis(~./125, name = "Mean Vote Share")) +
  ylab("Frequency") + xlab("All Candidates") +
  theme(axis.text.x = element_blank())
png("plot1_all.png", res = 75)
print(plot1.all)
dev.off()

plot1.gop <- ggplot(dat.gop.wfd, aes(X = wfd_ten)) +
  geom_col(aes(x = wfd_ten, y = gop.obs), size = 1, color = "gray65", fill = "white") +
  geom_line(aes(x = wfd_ten, y = 125*gop.mean), size = 1.5, color="gray65", group = 1) +
  scale_y_continuous(sec.axis = sec_axis(~./125, name = "Mean Vote Share")) +
  ylab("") + xlab("Republican Candidates") +
  theme(axis.text.x = element_blank())
png("plot1_gop.png", res = 75)
print(plot1.gop)
dev.off()

plot1.dem <- ggplot(dat.dem.wfd, aes(X = wfd_ten)) +
  geom_col(aes(x = wfd_ten, y = dem.obs), size = 1, color = "gray25", fill = "white") +
  geom_line(aes(x = wfd_ten, y = 125*dem.mean), size = 1.5, color="gray25", group = 1) +
  scale_y_continuous(sec.axis = sec_axis(~./125, name = "")) +
  ylab("Frequency") + xlab("Democratic Candidates") +
  theme(axis.text.x = element_blank())
png("plot1_dem.png", res = 75)
print(plot1.dem)
dev.off()

plot1.office.dem.out <- ggplot(dat.office.dem, aes(X = wfd_ten)) +
  geom_col(aes(x = wfd_ten, y = dem.obs), size = 1, color = "gray25", fill = "white") +
  geom_line(aes(x = wfd_ten, y = 125*dem.mean), size = 1.5, color = "gray25", group = 1) +
  scale_y_continuous(sec.axis = sec_axis(~./125, name = "")) +
  ylab("Frequency") + xlab("Democratic Candidates") +
  theme(axis.text.x = element_blank()) +
  facet_wrap(vars(OfficeCat))
png("plot1.office.dem.out.png", res = 75)
print(plot1.office.dem.out)
dev.off()

plot1.office.gop.out <- ggplot(dat.office.gop, aes(X = wfd_ten)) +
  geom_col(aes(x = wfd_ten, y = gop.obs), size = 1, color = "gray65", fill = "white") +
  geom_line(aes(x = wfd_ten, y = 125*gop.mean), size = 1.5, color = "gray65", group = 1) +
  scale_y_continuous(sec.axis = sec_axis(~./125, name = "")) +
  ylab("Frequency") + xlab("Republican Candidates") +
  theme(axis.text.x = element_blank()) +
  facet_wrap(vars(OfficeCat))
png("plot1.office.gop.out.png", res = 75)
print(plot1.office.gop.out)
dev.off()


## By race/party
plot1.dem.race <- ggplot(dat.race.dem, aes(X = wfd_ten)) +
  geom_col(aes(x = wfd_ten, y = dem.obs), size = 1, color = "gray20", fill = "white") +
  geom_line(aes(x = wfd_ten, y = 125*dem.mean), size = 1.5, color = "gray20", group = 1) +
  scale_y_continuous(sec.axis = sec_axis(~./125, name = "")) +
  ylab("") + xlab("Democratic Candidates") +
  theme(axis.text.x = element_blank()) +
  facet_wrap(vars(white))
png("plot1.dem.race.png", res = 75)
print(plot1.dem.race)
dev.off()

plot1.gop.race <- ggplot(dat.race.gop, aes(X = wfd_ten)) +
  geom_col(aes(x = wfd_ten, y = gop.obs), size = 1, color = "gray40", fill = "white") +
  geom_line(aes(x = wfd_ten, y = 125*gop.mean), size = 1.5, color = "gray40", group = 1) +
  scale_y_continuous(sec.axis = sec_axis(~./125, name = "")) +
  ylab("") + xlab("Republican Candidates") +
  theme(axis.text.x = element_blank()) +
  facet_wrap(vars(white))
png("plot1.gop.race.png", res = 75)
print(plot1.gop.race)
dev.off()


# Summary Table, Identity Variables and Outcome
table.identity <- data %>%
  dplyr::select(votepct, wfd_ten, CandidateDemocrat, white, black,
                latinohispanic, aapi, other, postgrad, governor, cabinet,
                otheroffice)
colnames(table.identity) <- c("Vote Pct", "Full WF Scale", "Cand. Dem", 
                              "Cand. White", "Cand. Black", 
                              "Cand. Latino", "Cand. AAPI", 
                              "Cand. Other Race", "Post-Grad Degree",
                              "Governor", "Cabinet Office", "Other Office")
mean.ti <- sapply(table.identity, mean, na.rm = T)
sd.ti <- sapply(table.identity, sd, na.rm = T)
obs.ti <- apply(table.identity, 2, function(x) sum((!is.na(x))))

tab.comb <- cbind(mean.ti, sd.ti, obs.ti)
colnames(tab.comb) <- c("Mean", "Std Dev", "Obs")


# Summary statistics, unique candidates
unique.table <- data %>%
  group_by(candid) %>%
  summarise(White = mean(white),
         Black = mean(black),
         Latino = mean(latinohispanic),
         AAPI = mean(aapi), OtherRace = mean(other),
         CandDem = mean(CandidateDemocrat)) %>%
  drop_na() %>% dplyr::select(-candid) %>%
  summarise_all(funs(sum))

party.table <- uniqcan %>%
  group_by(CandidateDemocrat, candid) %>%
  summarise(White = sum(white),
         Black = sum(black),
         Latino = sum(latinohispanic),
         AAPI = sum(aapi), OtherRace = sum(other)) %>%
  drop_na() %>% dplyr::select(-candid) %>%
  summarise_all(funs(sum))
```

```{r maps, echo = F, warning = F, message = F, include = F}
# Maps 
# Count of women by state
all_count <- data %>%
  group_by(State, State.FIPS) %>%
  summarise(count = n_distinct(candid)) %>%
  mutate(count = ifelse(State == "VA", NA, count))

# Manage shapefile cleaning
shapefile_counties <- get_urbn_map(map = "counties", sf = T)
shapefile_states <- get_urbn_map(map = "states", sf = T)

# Create map datasets for Dems, GOP #
dem_map <- datdem %>%
  group_by(fips) %>%
  summarise(vote_pct = mean(na.omit(votepct)))
gop_map <- datgop %>%
  group_by(fips) %>%
  summarise(vote_pct = mean(na.omit(votepct)))
all_map <- data %>%
  group_by(fips) %>%
  summarise(vote_pct = mean(na.omit(votepct)))
wfd3_map <- data %>%
  dplyr::select(fips, wfd_ten) %>%
  mutate(wfd_three = case_when(wfd_ten >= 0 & wfd_ten <= 2 ~ 0,
                               wfd_ten == 3 ~ 1,
                               wfd_ten >= 4 & wfd_ten <= 10 ~ 0,
                               .default = NA)) %>%
  group_by(fips) %>%
  summarise(wfd_three = mean(wfd_three, na.rm = T)) %>%
  mutate(wfd_three = case_when(wfd_three == 0 ~ "Not Three WFD Measures",
                               wfd_three == 1 ~ "Three WFD Measures",
                               .default = NA))


# Make maps
counts_map <- shapefile_states %>%
  left_join(., all_count, by = c("state_abbv" = "State")) %>%
  ggplot() + 
  geom_sf(aes(fill = count, color = "black", size = 0.25)) +
  geom_sf(data = shapefile_states, fill = NA, color = "black", size = 0.25) +
  scale_fill_gradient2(low = "gray95", mid = "gray55", high = "gray10",
                      na.value = "white", guide = "colourbar",
                      midpoint = 6.5,
                      breaks = c(0, 3, 6, 9, 12)) +
  scale_color_manual(values = c("black")) +
  theme(panel.background = element_blank(),
        axis.text.x = element_blank(), axis.text.y = element_blank(),
        axis.ticks = element_blank(), text = element_text(size=8),
        legend.position = "bottom",
        legend.key.size = unit(1.25, 'cm'),
        legend.text = element_text(size=10)) +
  labs(title = "Count of Female Candidates", fill = "") +
  ylab("") + xlab("") +
  guides(color = F, size = F, fill = guide_colorbar(nbin = 100))
png("counts_map.png", res = 75)
print(counts_map)
dev.off()

# county-wide maps
democrats_map <- shapefile_counties %>%
  mutate(county_fips = as.numeric(paste(county_fips, sep = ""))) %>%
  left_join(., dem_map, by = c("county_fips" = "fips")) %>%
  ggplot() +
  geom_sf(aes(fill = vote_pct, color = NA, size = 0.05)) +
  geom_sf(data = shapefile_states, fill = NA, color = "black", size = 0.25) +
  scale_fill_steps2(low = "gray95", mid = "gray45", high = "gray10",
                      na.value = "white", guide = "colourbar",
                      midpoint = 36.4,
                      breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) +
  scale_color_manual(values = c("black")) +
  theme(panel.background = element_blank(),
        axis.text.x = element_blank(), axis.text.y = element_blank(),
        axis.ticks = element_blank(), text = element_text(size=8),
        legend.position = "bottom",
        legend.key.size = unit(1.25, 'cm'),
        legend.text = element_text(size=10)) +
  labs(title = "Democratic Candidates", fill = "") +
  ylab("") + xlab("") +
  guides(color = F, size = F, fill = guide_colorbar(nbin = 100))
png("democrats_map.png", res = 75)
print(democrats_map)
dev.off()

republicans_map <- shapefile_counties %>%
  mutate(county_fips = as.numeric(paste(county_fips, sep = ""))) %>%
  left_join(., gop_map, by = c("county_fips" = "fips")) %>%
  ggplot() +
  geom_sf(aes(fill = vote_pct, color = NA, size = 0.05)) +
  geom_sf(data = shapefile_states, fill = NA, color = "black", size = 0.25) +
  scale_fill_steps2(low = "gray95", mid = "gray45", high = "gray10",
                      na.value = "white", guide = "colourbar",
                      midpoint = 36.4,
                      breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) +
  scale_color_manual(values = c("black")) +
  theme(panel.background = element_blank(),
        axis.text.x = element_blank(), axis.text.y = element_blank(),
        axis.ticks = element_blank(), text = element_text(size=8),
        legend.position = "bottom",
        legend.key.size = unit(1.25, 'cm'),
        legend.text = element_text(size=10)) +
  labs(title = "Republican Candidates", fill = "") +
  ylab("") + xlab("") +
  guides(color = F, size = F, fill = guide_colorbar(nbin = 100))
png("republicans_map.png", res = 75)
print(republicans_map)
dev.off()

candidates_map <- shapefile_counties %>%
  mutate(county_fips = as.numeric(paste(county_fips, sep = ""))) %>%
  left_join(., all_map, by = c("county_fips" = "fips")) %>%
  ggplot() +
  geom_sf(aes(fill = vote_pct, color = NA, size = 0.05)) +
  geom_sf(data = shapefile_states, fill = NA, color = "black", size = 0.25) +
  scale_fill_steps2(low = "gray95", mid = "gray45", high = "gray10",
                      na.value = "white", guide = "colourbar",
                      midpoint = 36.4,
                      breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) +
  scale_color_manual(values = c("black")) +
  theme(panel.background = element_blank(),
        axis.text.x = element_blank(), axis.text.y = element_blank(),
        axis.ticks = element_blank(), text = element_text(size=8),
        legend.position = "bottom",
        legend.key.size = unit(1.25, 'cm'),
        legend.text = element_text(size=10)) +
  labs(title = "All Candidates", fill = "") +
  ylab("") + xlab("") +
  guides(color = F, size = F, fill = guide_colorbar(nbin = 100))
png("candidates_map.png", res = 75)
print(candidates_map)
dev.off()

# Summary map, WFI == 3 only highlighted
wfdthree_map <- shapefile_counties %>%
  mutate(county_fips = as.numeric(paste(county_fips, sep = ""))) %>%
  left_join(., wfd3_map, by = c("county_fips" = "fips")) %>%
  ggplot() +
  geom_sf(aes(fill = wfd_three, color = NA, size = 0.05)) +
  geom_sf(data = shapefile_states, fill = NA, color = "black", size = 0.25) +
  scale_fill_manual(values = c("gray80", "gray20"), guide = "legend",
                      na.value = "white") +
  scale_color_manual(values = c("black")) +
  theme(panel.background = element_blank(),
        axis.text.x = element_blank(), axis.text.y = element_blank(),
        axis.ticks = element_blank(), text = element_text(size=8),
        legend.position = "bottom",
        legend.key.size = unit(1.25, 'cm'),
        legend.text = element_text(size=10)) +
  labs(title = "Counties with three measures of WFD index", fill = "") +
  ylab("") + xlab("") +
  guides(color = F, size = F, fill = guide_legend(ncol = 3))
png("wfdthree_map.png", res = 75)
print(wfdthree_map)
dev.off()

```

```{r models, echo = F, warning = F, message = F, include = F}
## Table, Aggregated Main Models with county-clustered SEs
ce1 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = data,
         se = "cluster", cluster = "fips")

ce2 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.demwh,
         se = "cluster", cluster = "fips")

ce3 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.demnw,
         se = "cluster", cluster = "fips")

ce4 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.gopwh,
         se = "cluster", cluster = "fips")

ce5 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.gopnw,
         se = "cluster", cluster = "fips")


# Testing out alternate metrics
## Table, eight components of WFD with county-clustered SEs
mw81 <- feols(votepct ~ manual_wfd_8 + wfd_8_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = data,
         se = "cluster", cluster = "fips")

mw82 <- feols(votepct ~ manual_wfd_8 + wfd_8_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.demwh,
         se = "cluster", cluster = "fips")

mw83 <- feols(votepct ~ manual_wfd_8 + wfd_8_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.demnw,
         se = "cluster", cluster = "fips")

mw84 <- feols(votepct ~ manual_wfd_8 + wfd_8_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.gopwh,
         se = "cluster", cluster = "fips")

mw85 <- feols(votepct ~ manual_wfd_8 + wfd_8_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.gopnw,
         se = "cluster", cluster = "fips")


## Table, eleven components of WFD with county-clustered SEs
mel1 <- feols(votepct ~ manual_wfd_11 + wfd_11_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = data,
         se = "cluster", cluster = "fips")

mel2 <- feols(votepct ~ manual_wfd_11 + wfd_11_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.demwh,
         se = "cluster", cluster = "fips")

mel3 <- feols(votepct ~ manual_wfd_11 + wfd_11_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.demnw,
         se = "cluster", cluster = "fips")

mel4 <- feols(votepct ~ manual_wfd_11 + wfd_11_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.gopwh,
         se = "cluster", cluster = "fips")

mel5 <- feols(votepct ~ manual_wfd_11 + wfd_11_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.gopnw,
         se = "cluster", cluster = "fips")

# Show distribution of the types of offices, broken by party and by office
governors_dem <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year), data = govdem,
         se = "cluster", cluster = "fips")

governors_gop <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year), data = govgop,
         se = "cluster", cluster = "fips")

cabinet_dem <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year), data = cabdem,
         se = "cluster", cluster = "fips")

cabinet_gop <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year), data = cabgop,
         se = "cluster", cluster = "fips")

other_dem <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year), data = othdem,
         se = "cluster", cluster = "fips")

other_gop <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year), data = othgop,
         se = "cluster", cluster = "fips")

# Analyse, democratic women of color, competitiveness
tr1 <- feols(votepct ~ wfd_ten + wfd_sq
          + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes, 
         data = dat.dwoc.safe,
         se = "cluster", cluster = "fips")

tr2 <- feols(votepct ~ wfd_ten + wfd_sq
          + victor + postgrad + mrp_ideology_mean
         + pct_dem + female_per
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes, 
         data = dat.dwoc.danger,
         se = "cluster", cluster = "fips")

tr3 <- feols(votepct ~ wfd_ten + wfd_sq
          + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes, 
         data = dat.dwoc.comp,
         se = "cluster", cluster = "fips")

# Analyse, democratic white women, competitiveness
tr4 <- feols(votepct ~ wfd_ten + wfd_sq
          + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes, 
         data = dat.dww.safe,
         se = "cluster", cluster = "fips")

tr5 <- feols(votepct ~ wfd_ten + wfd_sq
          + victor + postgrad + mrp_ideology_mean
         + pct_dem + female_per
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes, 
         data = dat.dww.danger,
         se = "cluster", cluster = "fips")

tr6 <- feols(votepct ~ wfd_ten + wfd_sq
          + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes, 
         data = dat.dww.comp,
         se = "cluster", cluster = "fips")

# Analyse, republican white women, competitiveness
tr7 <- feols(votepct ~ wfd_ten + wfd_sq
          + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes, 
         data = dat.rw.safe,
         se = "cluster", cluster = "fips")

tr8 <- feols(votepct ~ wfd_ten + wfd_sq
          + victor + postgrad + mrp_ideology_mean
         + pct_dem + female_per
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes, 
         data = dat.rw.danger,
         se = "cluster", cluster = "fips")

tr9 <- feols(votepct ~ wfd_ten + wfd_sq
          + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes, 
         data = dat.rw.comp,
         se = "cluster", cluster = "fips")
```

```{r matching, echo = F, warning = F, message = F, include = F}
## Match-as-preprocess, race

# Define vars
Y <- data$votepct
Tr <- data$white
X <- data %>% 
  dplyr::select(wfd_AAper, wfd_Hispper, wfd_foriegnb, wfd_medincome, 
                     wfd_college, wfd_bluecollar, wfd_schoolage, wfd_urban, 
                     wfd_nonsouth, wfd_age, CandidateDemocrat, victor, incumb, 
                     postgrad, mrp_ideology_mean, pct_dem, female_per, 
                     traditionalistic, moralistic, cathprop, evanprop, 
                     iwpr_rank, stateper, TotalVotes, Year, governor, cabinet)
ID <- data %>%
  dplyr::select(fips, OfficeCat)
t <- na.omit(cbind(Y, Tr, X, ID))
t_noY <- na.omit(cbind(Tr, X))

# Estimate balance on Tr
tr.predict <- glm(Tr ~ ., data = t_noY)
# summary(tr.predict)

Z <- t[,c(3:5, 7, 9, 11:16, 18, 20:25, 27:28)]

match <- Match(Y=t$Y, Tr=t$Tr, X=t[,c(3:29)], 
               Z=Z, BiasAdjust=T, estimand="ATT", 
               M=1, replace=F, caliper = 2.05)
# summary(match)

postmatch.dat <- as.data.frame(rbind(t[match$index.treated,],
                                     t[match$index.control,])) %>%
  as_tibble()
# tr.postdict <- glm(Tr ~ ., data = postmatch.dat[, -c(1)])
# summary(tr.postdict)

# Create new datasets
postmatch.dat$wfd_ten <- rowSums(postmatch.dat[,c(3:12)], na.rm = T)
postmatch.dat$wfd_sq <- I(postmatch.dat$wfd_ten^2)

match.whd <- postmatch.dat %>% 
  dplyr::filter(Tr == 1 & CandidateDemocrat == 1)
match.nwd <- postmatch.dat %>% 
  dplyr::filter(Tr == 0 & CandidateDemocrat == 1)
match.whr <- postmatch.dat %>% 
  dplyr::filter(Tr == 1 & CandidateDemocrat == 0)
match.nwr <- postmatch.dat %>% 
  dplyr::filter(Tr == 0 & CandidateDemocrat == 0)

## Matched dataset, models
mp1 <- feols(Y ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year), data = postmatch.dat,
         se = "cluster", cluster = "fips")

mp2 <- feols(Y ~ wfd_ten + wfd_sq
         + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = match.whd,
         se = "cluster", cluster = "fips")

mp3 <- feols(Y ~ wfd_ten + wfd_sq
         + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = match.nwd,
         se = "cluster", cluster = "fips")

mp4 <- feols(Y ~ wfd_ten + wfd_sq
         + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = match.whr,
         se = "cluster", cluster = "fips")

mp5 <- feols(Y ~ wfd_ten + wfd_sq
         + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = match.nwr,
         se = "cluster", cluster = "fips")
```

```{r marginaleffects1, echo = F, warning = F, message = F, include = F}
## Predicted effect models, unmatched data
wfd_uniq <- sort(unique(data$wfd_ten))
wfd_squniq <- sort(unique(data$wfd_sq))

# Make fitted marginal effects plot
WhiteDems <- ce2$coefficients[2] + (2 * ce2$coefficients[3] * wfd_squniq)
NonWhiteDems <- ce3$coefficients[2] + (2 * ce3$coefficients[3] * wfd_squniq)
WhiteGOPs <- ce4$coefficients[2] + (2 * ce4$coefficients[3] * wfd_squniq)
NonWhiteGOPs <- ce5$coefficients[2] + (2 * ce5$coefficients[3] * wfd_squniq)
wfd_plot <- as.data.frame(cbind(wfd_uniq,
                                WhiteDems, NonWhiteDems, 
                                WhiteGOPs, NonWhiteGOPs)) %>%
  gather(key = "variable", value = "value", -wfd_uniq)

# Create confidence intervals
iter <- 5000
sim.coefs.ce2 <- mvrnorm(n = iter, mu = ce2$coefficients, Sigma = vcov(ce2))
sim.coefs.ce3 <- mvrnorm(n = iter, mu = ce3$coefficients, Sigma = vcov(ce3))
sim.coefs.ce4 <- mvrnorm(n = iter, mu = ce4$coefficients, Sigma = vcov(ce4))
sim.coefs.ce5 <- mvrnorm(n = iter, mu = ce5$coefficients, Sigma = vcov(ce5))

# start by making marg effect for all models
marg.effect.data <- data.frame()
for(i in 1:iter){
  ME.ce2 <- sim.coefs.ce2[i,2] + 2*sim.coefs.ce2[i,3]*wfd_squniq
  ME.ce3 <- sim.coefs.ce3[i,2] + 2*sim.coefs.ce3[i,3]*wfd_squniq
  ME.ce4 <- sim.coefs.ce4[i,2] + 2*sim.coefs.ce4[i,3]*wfd_squniq
  ME.ce5 <- sim.coefs.ce5[i,2] + 2*sim.coefs.ce5[i,3]*wfd_squniq
  d <- data.frame(wfd_uniq, ME.ce2, ME.ce3, ME.ce4, ME.ce5)
  marg.effect.data <- bind_rows(marg.effect.data, d)
}

marg.effect.data <- marg.effect.data %>%
  group_by(wfd_uniq) %>%
  summarize(LB.whd = quantile(ME.ce2, .025),
            UB.whd = quantile(ME.ce2, .975),
            LB.nwd = quantile(ME.ce3, .025),
            UB.nwd = quantile(ME.ce3, .975),
            LB.whg = quantile(ME.ce4, .025),
            UB.whg = quantile(ME.ce4, .975),
            LB.nwg = quantile(ME.ce5, .025),
            UB.nwg = quantile(ME.ce5, .975)
            ) 

wfd_plot$LB <- as.numeric(paste(c(marg.effect.data$LB.whd, marg.effect.data$LB.nwd,
            marg.effect.data$LB.whg, marg.effect.data$LB.nwg), sep = ""))
wfd_plot$UB <- as.numeric(paste(c(marg.effect.data$UB.whd, marg.effect.data$UB.nwd,
            marg.effect.data$UB.whg, marg.effect.data$UB.nwg), sep = ""))

wfd_plot <- wfd_plot %>%
  mutate(sig = ifelse(sign(UB) == sign(LB), "significant", "not significant"))

# and plot
pd <- position_dodge(0.5)
plot.wfd.a <- ggplot(wfd_plot, aes(x = wfd_uniq, y = value, group = variable,
                                   color = sig)) +
  geom_errorbar(aes(ymin = LB, ymax = UB), width = 0.25, position = pd) +
  geom_point(position = pd) +
  geom_hline(yintercept = 0) +
  scale_color_manual(values = c("significant" = "gray10",
                                "not significant" = "gray60")) +
  labs(color = "", linetype = "", shape = "") +
  xlab("WFI Values") + ylab("Marginal Effect") + theme_bw() +
  theme(axis.text.x = element_blank(),
        legend.position = "none") +
  facet_wrap(facets = "variable", nrow = 2, ncol = 2)
```

```{r wfd_output1, message = F, warning = F, echo = F, include = F}
wfd_output <- rbind(tidy(governors_dem) %>%
                       dplyr::filter(term == "wfd_ten" | 
                                       term == "wfd_sq") %>%
                       mutate(model = "gov_dem"),
                    tidy(governors_gop) %>%
                       dplyr::filter(term == "wfd_ten" | 
                                       term == "wfd_sq") %>%
                       mutate(model = "gov_gop"),
                    tidy(cabinet_dem) %>%
                       dplyr::filter(term == "wfd_ten" | 
                                       term == "wfd_sq") %>%
                       mutate(model = "cab_dem"),
                    tidy(cabinet_gop) %>%
                       dplyr::filter(term == "wfd_ten" | 
                                       term == "wfd_sq") %>%
                       mutate(model = "cab_gop"),
                    tidy(other_dem) %>%
                       dplyr::filter(term == "wfd_ten" | 
                                       term == "wfd_sq") %>%
                       mutate(model = "oth_dem"),
                    tidy(other_gop) %>%
                       dplyr::filter(term == "wfd_ten" | 
                                       term == "wfd_sq") %>%
                       mutate(model = "oth_gop"))
```

```{r marginaleffects2, echo = F, warning = F, message = F, include = F}
## Predicted effect models, matched dataset models

# Create range of unique values
wfd_uniq <- sort(unique(data$wfd_ten))
wfd_squniq <- sort(unique(data$wfd_sq))

# Make fitted marginal effects plot
White_Dem <- mp2$coefficients[2] + (2 * mp2$coefficients[3] * wfd_squniq)
NonWhite_Dem <- mp3$coefficients[2] + (2 * mp3$coefficients[3] * wfd_squniq)
White_GOP <- mp4$coefficients[2] + (2 * mp4$coefficients[3] * wfd_squniq)
NonWhite_GOP <- mp5$coefficients[2] + (2 * mp5$coefficients[3] * wfd_squniq)

wfd_plot.mp <- as.data.frame(cbind(wfd_uniq,
                                White_Dem, NonWhite_Dem, 
                                White_GOP, NonWhite_GOP)) %>%
  gather(key = "variable", value = "value", -wfd_uniq)

# Create confidence intervals
iter <- 5000
sim.coefs.mp2 <- mvrnorm(n = iter, mu = mp2$coefficients, Sigma = vcov(mp2))
sim.coefs.mp3 <- mvrnorm(n = iter, mu = mp3$coefficients, Sigma = vcov(mp3))
sim.coefs.mp4 <- mvrnorm(n = iter, mu = mp4$coefficients, Sigma = vcov(mp4))
sim.coefs.mp5 <- mvrnorm(n = iter, mu = mp5$coefficients, Sigma = vcov(mp5))

# start by making marg effect for all models
marg.effect.match <- data.frame()
for(i in 1:iter){
  ME.mp2 <- sim.coefs.mp2[i,2] + 2*sim.coefs.mp2[i,3]*wfd_squniq
  ME.mp3 <- sim.coefs.mp3[i,2] + 2*sim.coefs.mp3[i,3]*wfd_squniq
  ME.mp4 <- sim.coefs.mp4[i,2] + 2*sim.coefs.mp4[i,3]*wfd_squniq
  ME.mp5 <- sim.coefs.mp5[i,2] + 2*sim.coefs.mp5[i,3]*wfd_squniq
  m <- data.frame(wfd_uniq, ME.mp2, ME.mp3, ME.mp4, ME.mp5)
  marg.effect.match <- bind_rows(marg.effect.match, m)
}

marg.effect.match <- marg.effect.match %>%
  group_by(wfd_uniq) %>%
  summarize(LB.whd = quantile(ME.mp2, .025),
            UB.whd = quantile(ME.mp2, .975),
            LB.nwd = quantile(ME.mp3, .025),
            UB.nwd = quantile(ME.mp3, .975),
            LB.whg = quantile(ME.mp4, .05),
            UB.whg = quantile(ME.mp4, .95),
            LB.nwg = quantile(ME.mp5, .05),
            UB.nwg = quantile(ME.mp5, .95)
            ) 

wfd_plot.mp$LB <- as.numeric(paste(c(marg.effect.match$LB.whd, marg.effect.match$LB.nwd,
            marg.effect.match$LB.whg, marg.effect.match$LB.nwg), sep = ""))
wfd_plot.mp$UB <- as.numeric(paste(c(marg.effect.match$UB.whd, marg.effect.match$UB.nwd,
            marg.effect.match$UB.whg, marg.effect.match$UB.nwg), sep = ""))

wfd_plot.mp <- wfd_plot.mp %>%
  mutate(sig = ifelse(sign(UB) == sign(LB), "significant", "not significant"))

# and plot
pd <- position_dodge(0.5)
plot.wfd.mp <- ggplot(wfd_plot.mp, aes(x = wfd_uniq, y = value, group = variable,
                                   color = sig)) +
  geom_errorbar(aes(ymin = LB, ymax = UB), width = 0.25, position = pd) +
  geom_point(position = pd) +
  geom_hline(yintercept = 0) +
  # geom_line(position = pd) +
  scale_color_manual(values = c("significant" = "gray10",
                                "not significant" = "gray60")) +
  labs(color = "", linetype = "", shape = "") +
  xlab("WFI Values") + ylab("Marginal Effect") + theme_bw() +
  theme(axis.text.x = element_blank(),
        legend.position = "none") +
  facet_wrap(facets = "variable", nrow = 2, ncol = 2)
```


```{r marginaleffects3, echo = F, warning = F, message = F, include = F}
## Predicted effect models, aggregated and matched models
wfd_uniq <- sort(unique(data$wfd_ten))
wfd_squniq <- sort(unique(data$wfd_sq))

govdem_out <- wfd_output$estimate[1] + (2 * wfd_output$estimate[2] * wfd_squniq)
govgop_out <- wfd_output$estimate[3] + (2 * wfd_output$estimate[4] * wfd_squniq)
cabdem_out <- wfd_output$estimate[5] + (2 * wfd_output$estimate[6] * wfd_squniq)
cabgop_out <- wfd_output$estimate[7] + (2 * wfd_output$estimate[8] * wfd_squniq)
othdem_out <- wfd_output$estimate[9] + (2 * wfd_output$estimate[10] * wfd_squniq)
othgop_out <- wfd_output$estimate[11] + (2 * wfd_output$estimate[12] * wfd_squniq)

wfd_plot <- as.data.frame(cbind(wfd_uniq, govdem_out, govgop_out, cabdem_out,
                                cabgop_out, othdem_out, othgop_out)) %>%
  gather(key = "variable", value = "value", -wfd_uniq) %>%
  mutate(variable = c(rep("Dem_Governor", 11), rep("GOP_Governor", 11),
                      rep("Dem_Cabinet", 11), rep("GOP_Cabinet", 11),
                      rep("Dem_Other", 11), rep("GOP_Other", 11)))

# Create confidence intervals
iter <- 5000
sim_govdem_out <- mvrnorm(n = iter, mu = governors_dem$coefficients, 
                          Sigma = vcov(governors_dem))
sim_govgop_out <- mvrnorm(n = iter, mu = governors_gop$coefficients, 
                          Sigma = vcov(governors_gop))
sim_cabdem_out <- mvrnorm(n = iter, mu = cabinet_dem$coefficients, 
                          Sigma = vcov(cabinet_dem))
sim_cabgop_out <- mvrnorm(n = iter, mu = cabinet_gop$coefficients, 
                          Sigma = vcov(cabinet_gop))
sim_othdem_out <- mvrnorm(n = iter, mu = other_dem$coefficients, 
                          Sigma = vcov(other_dem))
sim_othgop_out <- mvrnorm(n = iter, mu = other_gop$coefficients, 
                          Sigma = vcov(other_gop))

marg.effect.data <- data.frame()
for(i in 1:iter){
  ME.govdem <- sim_govdem_out[i,2] + 2*sim_govdem_out[i,3]*wfd_squniq
  ME.govgop <- sim_govgop_out[i,2] + 2*sim_govgop_out[i,3]*wfd_squniq
  ME.cabdem <- sim_cabdem_out[i,2] + 2*sim_cabdem_out[i,3]*wfd_squniq
  ME.cabgop <- sim_cabgop_out[i,2] + 2*sim_cabgop_out[i,3]*wfd_squniq
  ME.othdem <- sim_othdem_out[i,2] + 2*sim_othdem_out[i,3]*wfd_squniq
  ME.othgop <- sim_othgop_out[i,2] + 2*sim_othgop_out[i,3]*wfd_squniq
  d <- data.frame(wfd_uniq, ME.govdem, ME.govgop, ME.cabdem, ME.cabgop,
                  ME.othdem, ME.othgop)
  marg.effect.data <- bind_rows(marg.effect.data, d)
}

marg.effect.data <- marg.effect.data %>%
  group_by(wfd_uniq) %>%
  summarise(lb.govdem = quantile(ME.govdem, 0.025),
            ub.govdem = quantile(ME.govdem, 0.975),
            lb.govgop = quantile(ME.govgop, 0.025),
            ub.govgop = quantile(ME.govgop, 0.975),
            lb.cabdem = quantile(ME.cabdem, 0.025),
            ub.cabdem = quantile(ME.cabdem, 0.975),
            lb.cabgop = quantile(ME.cabgop, 0.025),
            ub.cabgop = quantile(ME.cabgop, 0.975),
            lb.othdem = quantile(ME.othdem, 0.025),
            ub.othdem = quantile(ME.othdem, 0.975),
            lb.othgop = quantile(ME.othgop, 0.025),
            ub.othgop = quantile(ME.othgop, 0.975)
            )

wfd_plot$LB <- as.numeric(paste(c(marg.effect.data$lb.govdem, marg.effect.data$lb.govgop,
                                  marg.effect.data$lb.cabdem, marg.effect.data$lb.cabgop,
                                  marg.effect.data$lb.othdem, marg.effect.data$lb.othgop)))
wfd_plot$UB <- as.numeric(paste(c(marg.effect.data$ub.govdem, marg.effect.data$ub.govgop,
                                  marg.effect.data$ub.cabdem, marg.effect.data$ub.cabgop,
                                  marg.effect.data$ub.othdem, marg.effect.data$ub.othgop)))
wfd_plot$party <- c(rep(c(rep("Dem", 11), rep("GOP", 11)), 3))

wfd_plot <- wfd_plot %>%
  mutate(sig = ifelse(sign(UB) == sign(LB), "significant", "not significant"))

# and plot
pd <- position_dodge(0.75)
override.shape = c(16, 16, 18, 18, 20, 20)
override.linetype = c(1, 3, 1, 3, 1, 3)

plot.wfd.office <- ggplot(wfd_plot, aes(x = wfd_uniq, y = value, 
                                    color = sig,
                                    group = variable,
                                    shape = party,
                                    linetype = variable)) +
  geom_errorbar(aes(ymin = LB, ymax = UB), width = 0.25, position = pd) +
  geom_point(position = pd) +
  geom_hline(yintercept = 0) +
  labs(color = "", linetype = "", shape = "") +
  xlab("WFI Values") + ylab("Marginal Effect") + theme_bw() +
  theme(axis.text.x = element_blank(), legend.position = "none") +
  # guides(color = guide_legend(override.aes = list(shape = override.shape, 
  #                                                 linetype = override.linetype))) +
  scale_shape(guide = "none") + 
  scale_linetype(guide = "none") +
  scale_color_manual(values = c("significant" = "gray10",
                                "not significant" = "gray60")) +
  facet_wrap(facets = "variable", ncol = 3, nrow = 2)
```

```{r marginaleffects4, echo = F, warning = F, message = F, include = F}
## Predicted effect models, truncated Dem women of color models
wfd_uniq <- sort(unique(data$wfd_ten))
wfd_squniq <- sort(unique(data$wfd_sq))

# Make fitted marginal effects plot, truncated dataset
SafeDem <- tr1$coefficients[2] + (2 * tr1$coefficients[3] * wfd_squniq)
SafeGOP <- tr2$coefficients[2] + (2 * tr2$coefficients[3] * wfd_squniq)
Competitive <- tr3$coefficients[2] + (2 * tr3$coefficients[3] * wfd_squniq)
trn_plot <- as.data.frame(cbind(wfd_uniq,
                                SafeDem, SafeGOP,
                                Competitive)) %>%
  gather(key = "variable", value = "value", -wfd_uniq)

# Create confidence intervals
iter <- 5000
sim.coefs.tr1 <- mvrnorm(n = iter, mu = tr1$coefficients, Sigma = vcov(tr1))
sim.coefs.tr2 <- mvrnorm(n = iter, mu = tr2$coefficients, Sigma = vcov(tr2))
sim.coefs.tr3 <- mvrnorm(n = iter, mu = tr3$coefficients, Sigma = vcov(tr3))

# start by making marg effect for all models
marg.effect.data <- data.frame()
for(i in 1:iter){
  ME.tr1 <- sim.coefs.tr1[i,2] + 2*sim.coefs.tr1[i,3]*wfd_squniq
  ME.tr2 <- sim.coefs.tr2[i,2] + 2*sim.coefs.tr2[i,3]*wfd_squniq
  ME.tr3 <- sim.coefs.tr3[i,2] + 2*sim.coefs.tr3[i,3]*wfd_squniq
  d <- data.frame(wfd_uniq, ME.tr1, ME.tr2, ME.tr3)
  marg.effect.data <- bind_rows(marg.effect.data, d)
}

marg.effect.data <- marg.effect.data %>%
  group_by(wfd_uniq) %>%
  summarize(LB.tr1 = quantile(ME.tr1, .05),
            UB.tr1 = quantile(ME.tr1, .95),
            LB.tr2 = quantile(ME.tr2, .05),
            UB.tr2 = quantile(ME.tr2, .95),
            LB.tr3 = quantile(ME.tr3, .05),
            UB.tr3 = quantile(ME.tr3, .95)
            ) 

trn_plot$LB <- as.numeric(paste(c(marg.effect.data$LB.tr1, marg.effect.data$LB.tr2,
            marg.effect.data$LB.tr3), sep = ""))
trn_plot$UB <- as.numeric(paste(c(marg.effect.data$UB.tr1, marg.effect.data$UB.tr2,
            marg.effect.data$UB.tr3), sep = ""))

trn_plot <- trn_plot %>%
  mutate(sig = ifelse(sign(UB) == sign(LB), "significant", "not significant"))

# and plot
pd <- position_dodge(0.5)
plot.wfd.tr <- ggplot(trn_plot, aes(x = wfd_uniq, y = value, group = variable,
                                   color = sig)) +
  geom_errorbar(aes(ymin = LB, ymax = UB), width = 0.25, position = pd) +
  geom_point(position = pd) +
  geom_hline(yintercept = 0) +
  # geom_line(position = pd) +
  scale_color_manual(values = c("significant" = "gray10",
                                "not significant" = "gray60")) +
  labs(color = "", linetype = "", shape = "") +
  xlab("WFI Values") + ylab("Marginal Effect") + theme_bw() +
  theme(axis.text.x = element_blank(),
        legend.position = "none") +
  facet_wrap(facets = "variable", nrow = 1, ncol = 3)
```

```{r marginaleffects5, echo = F, warning = F, message = F, include = F}
## Predicted effect models, truncated Dem white women models
wfd_uniq <- sort(unique(data$wfd_ten))
wfd_squniq <- sort(unique(data$wfd_sq))

# Make fitted marginal effects plot, truncated dataset
SafeDem <- tr4$coefficients[2] + (2 * tr4$coefficients[3] * wfd_squniq)
SafeGOP <- tr5$coefficients[2] + (2 * tr5$coefficients[3] * wfd_squniq)
Competitive <- tr6$coefficients[2] + (2 * tr6$coefficients[3] * wfd_squniq)
trn_plot <- as.data.frame(cbind(wfd_uniq,
                                SafeDem, SafeGOP,
                                Competitive)) %>%
  gather(key = "variable", value = "value", -wfd_uniq)

# Create confidence intervals
iter <- 5000
sim.coefs.tr4 <- mvrnorm(n = iter, mu = tr4$coefficients, Sigma = vcov(tr4))
sim.coefs.tr5 <- mvrnorm(n = iter, mu = tr5$coefficients, Sigma = vcov(tr5))
sim.coefs.tr6 <- mvrnorm(n = iter, mu = tr6$coefficients, Sigma = vcov(tr6))

# start by making marg effect for all models
marg.effect.data <- data.frame()
for(i in 1:iter){
  ME.tr4 <- sim.coefs.tr4[i,2] + 2*sim.coefs.tr4[i,3]*wfd_squniq
  ME.tr5 <- sim.coefs.tr5[i,2] + 2*sim.coefs.tr5[i,3]*wfd_squniq
  ME.tr6 <- sim.coefs.tr6[i,2] + 2*sim.coefs.tr6[i,3]*wfd_squniq
  d <- data.frame(wfd_uniq, ME.tr4, ME.tr5, ME.tr6)
  marg.effect.data <- bind_rows(marg.effect.data, d)
}

marg.effect.data <- marg.effect.data %>%
  group_by(wfd_uniq) %>%
  summarize(LB.tr4 = quantile(ME.tr4, .05),
            UB.tr4 = quantile(ME.tr4, .95),
            LB.tr5 = quantile(ME.tr5, .05),
            UB.tr5 = quantile(ME.tr5, .95),
            LB.tr6 = quantile(ME.tr6, .05),
            UB.tr6 = quantile(ME.tr6, .95)
            ) 

trn_plot$LB <- as.numeric(paste(c(marg.effect.data$LB.tr4, marg.effect.data$LB.tr5,
            marg.effect.data$LB.tr6), sep = ""))
trn_plot$UB <- as.numeric(paste(c(marg.effect.data$UB.tr4, marg.effect.data$UB.tr5,
            marg.effect.data$UB.tr6), sep = ""))

trn_plotw <- trn_plot %>%
  mutate(sig = ifelse(sign(UB) == sign(LB), "significant", "not significant"))

# and plot
pd <- position_dodge(0.5)
plot.wfd.trw <- ggplot(trn_plotw, aes(x = wfd_uniq, y = value, group = variable,
                                   color = sig)) +
  geom_errorbar(aes(ymin = LB, ymax = UB), width = 0.25, position = pd) +
  geom_point(position = pd) +
  geom_hline(yintercept = 0) +
  # geom_line(position = pd) +
  scale_color_manual(values = c("significant" = "gray10",
                                "not significant" = "gray60")) +
  labs(color = "", linetype = "", shape = "") +
  xlab("WFI Values") + ylab("Marginal Effect") + theme_bw() +
  theme(axis.text.x = element_blank(),
        legend.position = "none") +
  facet_wrap(facets = "variable", nrow = 1, ncol = 3)
```

```{r marginaleffects6, echo = F, warning = F, message = F, include = F}
## Predicted effect models, truncated Republican women models
wfd_uniq <- sort(unique(data$wfd_ten))
wfd_squniq <- sort(unique(data$wfd_sq))

# Make fitted marginal effects plot, truncated dataset
SafeDem <- tr7$coefficients[2] + (2 * tr7$coefficients[3] * wfd_squniq)
SafeGOP <- tr8$coefficients[2] + (2 * tr8$coefficients[3] * wfd_squniq)
Competitive <- tr9$coefficients[2] + (2 * tr9$coefficients[3] * wfd_squniq)
trn_plot <- as.data.frame(cbind(wfd_uniq,
                                SafeDem, SafeGOP,
                                Competitive)) %>%
  gather(key = "variable", value = "value", -wfd_uniq)

# Create confidence intervals
iter <- 5000
sim.coefs.tr7 <- mvrnorm(n = iter, mu = tr7$coefficients, Sigma = vcov(tr7))
sim.coefs.tr8 <- mvrnorm(n = iter, mu = tr8$coefficients, Sigma = vcov(tr8))
sim.coefs.tr9 <- mvrnorm(n = iter, mu = tr9$coefficients, Sigma = vcov(tr9))

# start by making marg effect for all models
marg.effect.data <- data.frame()
for(i in 1:iter){
  ME.tr7 <- sim.coefs.tr7[i,2] + 2*sim.coefs.tr7[i,3]*wfd_squniq
  ME.tr8 <- sim.coefs.tr8[i,2] + 2*sim.coefs.tr8[i,3]*wfd_squniq
  ME.tr9 <- sim.coefs.tr9[i,2] + 2*sim.coefs.tr9[i,3]*wfd_squniq
  d <- data.frame(wfd_uniq, ME.tr7, ME.tr8, ME.tr9)
  marg.effect.data <- bind_rows(marg.effect.data, d)
}

marg.effect.data <- marg.effect.data %>%
  group_by(wfd_uniq) %>%
  summarize(LB.tr7 = quantile(ME.tr7, .05),
            UB.tr7 = quantile(ME.tr7, .95),
            LB.tr8 = quantile(ME.tr8, .05),
            UB.tr8 = quantile(ME.tr8, .95),
            LB.tr9 = quantile(ME.tr9, .05),
            UB.tr9 = quantile(ME.tr9, .95)
            ) 

trn_plot$LB <- as.numeric(paste(c(marg.effect.data$LB.tr7, marg.effect.data$LB.tr8,
            marg.effect.data$LB.tr9), sep = ""))
trn_plot$UB <- as.numeric(paste(c(marg.effect.data$UB.tr7, marg.effect.data$UB.tr8,
            marg.effect.data$UB.tr9), sep = ""))

trn_plotr <- trn_plot %>%
  mutate(sig = ifelse(sign(UB) == sign(LB), "significant", "not significant"))

# and plot
pd <- position_dodge(0.5)
plot.wfd.trr <- ggplot(trn_plotr, aes(x = wfd_uniq, y = value, group = variable,
                                   color = sig)) +
  geom_errorbar(aes(ymin = LB, ymax = UB), width = 0.25, position = pd) +
  geom_point(position = pd) +
  geom_hline(yintercept = 0) +
  # geom_line(position = pd) +
  scale_color_manual(values = c("significant" = "gray10",
                                "not significant" = "gray60")) +
  labs(color = "", linetype = "", shape = "") +
  xlab("WFI Values") + ylab("Marginal Effect") + theme_bw() +
  theme(axis.text.x = element_blank(),
        legend.position = "none") +
  facet_wrap(facets = "variable", nrow = 1, ncol = 3)
```


```{r margplot2_setup, include=FALSE, echo=FALSE}
varx <- seq(from = 1, to = 10, by = 1)
vary1 <- seq(from = 0.1, to = 2.0, by = 0.2)
vary2 <- seq(from = 0.5, to = 1.4, by = 0.1)
vary3 <- seq(from = 2.1, to = 4.0, by = 0.2)
vary4 <- seq(from = 2.5, to = 3.4, by = 0.1)
```

```{r margplot2, echo=FALSE, fig.cap="Corresponds to Figure 1 in manuscript", out.width = "100%"}
plot(vary1 ~ varx, type = "l", lwd = 3, lty = 1,
     ylim = c(0, 6.1), col = "gray55",
     xlab = "Women-Friendliness", ylab = "Electoral Performance",
     xaxt = "n", yaxt = "n", main = "")
lines(vary2 ~ varx, type = "l", lwd = 3, lty = 2, col = "gray55")
lines(vary3 ~ varx, type = "l", lwd = 3, lty = 1, col = "gray5")
lines(vary4 ~ varx, type = "l", lwd = 3, lty = 2, col = "gray5")
legend(0.35, 6.3, legend = c("Dem (Non-White)", "Dem (White)", 
                             "GOP (Non-White)", "GOP (White)"),
       col = c("gray5", "gray5", "gray55", "gray55"),
       lty = c(1,2,1,2), cex = 0.95, lwd = 3)
```

```{r candidates_count, fig.show="hold", echo=FALSE, warning=FALSE, fig.cap="Corresponds to Figure 2 in manuscript", out.width = "100%"}
plot(counts_map)
```

```{r candidates_map, fig.show="hold", echo=FALSE, warning=FALSE, fig.cap="Corresponds to Figure 3 in manuscript", out.width = "100%"}
plot(candidates_map)
```

```{r summary_race_info, echo = F, warning = F, message = F, results = "asis", comment = "", fig.cap = "Corresponds to Table 1 in manuscript"}
race_stats <- data %>%
  group_by(candid) %>%
  summarise(white = mean(white, na.rm = T),
            dem = mean(CandidateDemocrat, na.rm = T)) %>%
  ungroup() %>%
  mutate(dem = round(dem, 0),
         dem = case_when(dem == 1 ~ "Democrat", .default = "Republican"),
         white = case_when(white == 1 ~ "White", .default = "Non-White"))
options(width = 60)
table_ss <- table(race_stats$white, race_stats$dem)
kable(table_ss, caption = "",
      digits = 3)
```

```{r plot1.party1, figures-side, fig.show="hold", out.width="100%", echo=FALSE, warning=FALSE, fig.cap="Corresponds to Figure 4 in manuscript"}
grid.arrange(plot1.dem, plot1.gop, ncol = 1, nrow = 2)
```

```{r plot1.party2, figures-side, fig.show="hold", out.width="100%", echo=FALSE, warning=FALSE, fig.cap="Corresponds to Figure 5 in manuscript"}
grid.arrange(plot1.dem.race, plot1.gop.race, ncol = 1, nrow = 2)
```

```{r modelsummary_1, comment = "", echo = F, message = F, warning = F, fig.cap = "Corresponds to Table 2 in manuscript"}
models_ce1_ce5 <- list("All" = ce1, "White Democrats" = ce2, "NonWhite Democrats" = ce3, 
          "White Republicans" = ce4, "NonWhite Republicans" = ce5)
rsquares_ce1_ce5 <- c("Adj. r2", 
                      round(as.numeric(fitstat(ce1, "ar2")), 3), 
                      round(as.numeric(fitstat(ce2, "ar2")), 3), 
                      round(as.numeric(fitstat(ce3, "ar2")), 3), 
                      round(as.numeric(fitstat(ce4, "ar2")), 3), 
                      round(as.numeric(fitstat(ce5, "ar2")), 3))
rmse_ce1_ce5 <- c("RMSE", 
                      round(as.numeric(fitstat(ce1, "rmse")), 2), 
                      round(as.numeric(fitstat(ce2, "rmse")), 2), 
                      round(as.numeric(fitstat(ce3, "rmse")), 2), 
                      round(as.numeric(fitstat(ce4, "rmse")), 2), 
                      round(as.numeric(fitstat(ce5, "rmse")), 2))
rows_ce1_ce5 <- as.data.frame(rbind(rsquares_ce1_ce5, rmse_ce1_ce5))

modelsummary(models_ce1_ce5, output = "markdown",
             fmt_statistic(estimate = 3),
             coef_omit = c(-1, -2, -3),
             vcovBS = "bootstrap", R = 1000, stars = T,
             coef_rename = c("(Intercept)" = "Intercept", "wfd_ten" = "WFI (10)",
                             "wfd_sq" = "WFI Squared"),
             gof_map = c("nobs", "r2.adjusted", "F"), align = "lddddd",
             add_rows = rows_ce1_ce5,
             title = "The Effect of Women-Friendliness on U.S. Statewide Executive Woman Candidate County-Level Vote Share, 2010-2019. Robust standard errors reported in parentheses.")
```

```{r plot_wfd_a, echo=FALSE, warning=FALSE, fig.cap="Corresponds to Figure 6 in manuscript", fig.show="hold", out.width="100%"}
plot.wfd.a
```

```{r modelsummary_2, comment = "", echo = F, message = F, warning = F, fig.cap = "Corresponds to Table 3 in manuscript"}
models_mp1_mp5 <- list("All" = mp1, "White Democrats" = mp2, "NonWhite Democrats" = mp3, 
          "White Republicans" = mp4, "NonWhite Republicans" = mp5)
rsquares_mp1_mp5  <- c("Adj. r2", 
                      round(as.numeric(fitstat(mp1, "ar2")), 3), 
                      round(as.numeric(fitstat(mp2, "ar2")), 3), 
                      round(as.numeric(fitstat(mp3, "ar2")), 3), 
                      round(as.numeric(fitstat(mp4, "ar2")), 3), 
                      round(as.numeric(fitstat(mp5, "ar2")), 3))
rmse_mp1_mp5  <- c("RMSE", 
                      round(as.numeric(fitstat(mp1, "rmse")), 2), 
                      round(as.numeric(fitstat(mp2, "rmse")), 2), 
                      round(as.numeric(fitstat(mp3, "rmse")), 2), 
                      round(as.numeric(fitstat(mp4, "rmse")), 2), 
                      round(as.numeric(fitstat(mp5, "rmse")), 2))
rows_mp1_mp5  <- as.data.frame(rbind(rsquares_mp1_mp5, rmse_mp1_mp5))

modelsummary(models_mp1_mp5, output = "markdown",
             fmt_statistic(estimate = 3),
             coef_omit = c(-1, -2, -3),
             vcovBS = "bootstrap", R = 1000, stars = T,
             coef_rename = c("(Intercept)" = "Intercept", "wfd_ten" = "WFI (10)",
                             "wfd_sq" = "WFI Squared"),
             gof_map = c("nobs", "r2.adjusted", "F"), align = "lddddd",
             add_rows = rows_mp1_mp5, 
             title = "The Effect of Women-Friendliness on U.S. Statewide Executive Woman Candidate County-Level Vote Share, 2010-2019. Robust standard errors reported in parentheses.")
```

```{r plot_wfd_mp, echo=FALSE, warning=FALSE, fig.cap="Corresponds to Figure 7 in manuscript", fig.show="hold", out.width="100%"}
plot.wfd.mp
```

```{r plot_wfd_tr, echo=FALSE, warning=FALSE, fig.cap="Corresponds to Figure 8 in manuscript"}
plot.wfd.tr
plot.wfd.trw
plot.wfd.trr
```


